%% Example of Single Molecule Localization


%Make an image of a single molecule observed by a microscope

Pixels=8;
Sigma=1.5;
Photons = 400;
Background = 1;

[XGrid,YGrid]=meshgrid((1:Pixels),(1:Pixels))

%Use Gaussian approximation to the Point Spread Function
Model=@(Theta)Theta(3)*normpdf(XGrid,Theta(2),Sigma).*normpdf(YGrid,Theta(1),Sigma)+Theta(4);


X=4;
Y=5;
TrueData = Model([Y,X,Photons,Background])
figure
imshow(TrueData,[],'InitialMagnification','fit')

Data = poissrnd(TrueData)
figure
imshow(Data,[],'InitialMagnification','fit')

%% Markov-Chain Monte-Carlo

N=10000;  %number of steps in  chain
 
% Parameters are [Y X Photons BG];
%Intitialize chain vector
Chain=zeros(N,4);

%Setup initial value
Chain(1,:)=[4, 4, max(Data(:)), min(Data(:))];

%Setup jump sizes for each parameter
SigmaJump=[.1 .1 10 .1];

%run chain

for nn=1:N-1
    %Make test jump
    CurrentX=Chain(nn,:);
    TestX=CurrentX+SigmaJump.*randn(1,4);
    
    %Calculate the likelihood ratio
    ModelCurrent = Model(CurrentX);
    ModelTest = Model(TestX);
    
    %Log Likelihood Ratio
    LLR = - ModelTest + ModelCurrent + Data.*(log(ModelTest)-log(ModelCurrent));
    LLR=sum(LLR(:));
    LikelihoodRatio = exp(LLR);
    
    %Metropolis-Hastings rule
    if rand()< LikelihoodRatio
        Chain(nn+1,:)=TestX;
    else
        Chain(nn+1,:)=CurrentX;
    end
end

figure
plot(Chain(:,2),Chain(:,1))
xlabel('X (Pixels)')
ylabel('Y (Pixels)')


figure
plot(Chain(:,2),Chain(:,1),'.')
xlabel('X (Pixels)')
ylabel('Y (Pixels)')

%Find mean and variance
NBurnin = 1000;
ThetaMean = mean(Chain(NBurnin:end,:))
ThetaSE = std(Chain(NBurnin:end,:))


